package Math::Geometry::ConeMap;
$VERSION = v0.0.1;

use warnings;
use strict;
use Carp;

use utf8;

sub π () {atan2(1,1) * 4};

=encoding utf8

=head1 NAME

Math::Geometry::ConeMap - mapping 2D points onto a cone

=head1 SYNOPSIS

The cone geometry is described by the apex angle and the mapping is
governed by a distance from the apex.

The apex angle (C<θ>) is given in radians (this is the complement of the
base angle.)

The distance (C<D>) is measured from the apex.  This describes a
waterline on the cone at which the circumference corresponds to linear
measurement on the x-axis of the input data.  That is, the input points
will be distorted (compressed) above and (stretched) below this
waterline.

  my $cone = Math::Geometry::ConeMap->new(θ => π/4, D => 3.2);

  my $fD = $cone->fD_of_y;
  my $fφ = $cone->fφ_of_x;
  my $fxyz1 = $cone->xyz_of_xy;
  my $fxyz2 = $cone->xyz_of_Dφ;

=head1 Bounds

The input y-values should not exceed C<D>.  Values above this correspond
to points above the cone apex.  This implementation may not support such
imaginary "mirrored cone" surfaces.

The absolute value of input x-values should not exceed C<π*D*sin(θ)>.
Values outside of this range will result in C<φ> values which are
greater than C<π> away from the baseline (i.e. wound around the cone
more than once.)

=cut

=head1 Constructor

=head2 new

  my $cone = Math::Geometry::ConeMap->new(θ => π/4, D => 3.2);

=cut

sub new {
  my $class = shift;
  my $self = {@_};
  bless($self, $class);

  croak("must have a centerline distance") unless(exists $self->{D});
  croak("must have an apex angle")         unless(exists $self->{θ});

  return($self);
} # end subroutine new definition
########################################################################

sub _function {
  my ($of, $str) = @_;

  $str =~ s/(\w+?)²/($1*$1)/g; # squaring

  foreach my $rep ($str =~ m/^\s*my (\w+) =/mg) {
    #warn "replace $rep";
    $str =~ s/\b$rep\b/\$$rep/g;
  }

  my @of = split(//, $of);
  foreach my $n (0..$#of) {
    $str =~ s/\b$of[$n]\b/\$_[$n]/g;
  }
  my $sub = 'sub ('. '$'x@of . ') {'. $str . '}';
  # warn "now compile <<$sub>>\n";
  _compile($sub);
}
sub _compile {
  my $sub = eval($_[0]);
  $@ and croak("compiling '$_[0]' failed: $@");
  return($sub);
}

=head1 Functions

The cone parameters D and θ are constant for any given cone object.
Therefore, speed says we have function generators here instead of
methods.  Each of the following methods will return a compiled function.

=head2 D_of_y

Returns a function which calculates the distance D from the cone's apex
given a y-value as input.

  my $fD = $cone->D_of_y;

  my $D = $fD->($y);

=cut

sub D_of_y {
  my $self = shift;

  my $D = $self->{D};
  return _function(Y => "$D - Y");
} # end subroutine D_of_y definition
########################################################################

=head2 φ_of_x

Returns a function which calculates the angle φ about the cone's axis
given an x-value as input.  Zero is taken as the x-axis with positive
angles counter-clockwise (i.e. the z-axis corresponds to the cone axis.)

  my $fφ = $cone->φ_of_x;

  my $φ = $fφ->($x);

=cut

sub φ_of_x {
  my $self = shift;

  my $D = $self->{D};
  my $θ = $self->{θ};
  my $div = $D * sin($θ);

  return _function(X => "X/$div");
} # end subroutine φ_of_x definition
########################################################################

=head2 xyz_of_xy

Returns a function which maps a 3D point on the cone surface from an
x-value,y-value input.

  my $fxyz = $cone->xyz_of_xy;

  my @xy = $fxyz->($x, $y);

=cut

sub xyz_of_xy {
  my $self = shift;

  my $D = $self->{D};
  my $θ = $self->{θ};
  my $div = $D * sin($θ);

  my $sinθ = sin($θ);
  my $cosθ = cos($θ);

  return _function(XY => "
    my d = $D - Y; # fD(y)
    my φ = X/$div; # fφ(x)
    my rxy = d * $sinθ;
    return(
      rxy * cos(φ), # x
      rxy * sin(φ), # y
      -d  * $cosθ,  # z
  )");
} # end subroutine xyz_of_xy definition
########################################################################

=head2 xyz_of_Dφ

Returns a function which calculates a 3D point on the cone surface given
a distance D and angle φ.

  my $fxyz = $cone->xyz_of_Dφ;

  my ($x, $y, $z) = $fxyz->($D, $φ);

=cut

sub xyz_of_Dφ {
  my $self = shift;

  my $Dcl = $self->{D};
  my $θ = $self->{θ};
  my $sinθ = sin($θ);
  my $cosθ = cos($θ);

  return _function(dφ => "
    my rxy = d * $sinθ;
    return(
      rxy * cos(φ), # x
      rxy * sin(φ), # y
      -d * $cosθ,   # z
    )");
} # end subroutine xyz_of_Dφ definition
########################################################################

=head2 z_of_R

Calculate the z-value of a point on the cone at perpendicular distance
C<$R> from the axis.

  my $fZ = $cone->z_of_R();

  my $z = $fZ->($R);
=cut

sub z_of_R {
  my $self = shift;

  my $θ = $self->{θ};
  my $tanθ = sin($θ) / cos($θ);
  return _function(R => "-R/$tanθ");
} # end subroutine z_of_R definition
########################################################################

=head2 z_of_xy

Same as z_of_R(), but includes ƒR(XY) = sqrt(X²+Y²).  The X,Y here is
within 3D space.

  my $fZ = $cone->z_of_xy();

  my $z = $fZ->($X, $Y);

=cut

sub z_of_xy {
  my $self = shift;

  my $θ = $self->{θ};
  my $tanθ = sin($θ) / cos($θ);
  return _function(XY => "-sqrt(X²+Y²)/$tanθ");
} # end subroutine z_of_xy definition
########################################################################

=head1 AUTHOR

Eric Wilhelm @ <ewilhelm at cpan dot org>

http://scratchcomputing.com/

=head1 BUGS

If you found this module on CPAN, please report any bugs or feature
requests through the web interface at L<http://rt.cpan.org>.  I will be
notified, and then you'll automatically be notified of progress on your
bug as I make changes.

If you pulled this development version from my /svn/, please contact me
directly.

=head1 COPYRIGHT

Copyright (C) 2009 Eric L. Wilhelm, All Rights Reserved.

=head1 NO WARRANTY

Absolutely, positively NO WARRANTY, neither express or implied, is
offered with this software.  You use this software at your own risk.  In
case of loss, no person or entity owes you anything whatsoever.  You
have been warned.

=head1 LICENSE

This program is free software; you can redistribute it and/or modify it
under the same terms as Perl itself.

=cut

# vi:ts=2:sw=2:et:sta
1;
